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Abstract 


A method for using integrating detectors to estimate 
the phase of an optical wavefront is investigated. The 
phase and amplitude are assumed to vary slowly compared 
to the integration time and the integrated samples are 
shown to be corrupted by white gaussian noise. A maxi¬ 
mum a-posteriori nonlinear Kalman filter is derived and 
simulated for both constant and random walk phase pro¬ 
cesses. The performance of the filter is compared to its 
Cramer-Rao lower bound and a first order linearized phase- 
locked loop (PLL). The filter never performs a great deal 
better than the PLL, but it can be implemented in a low 
bandwidth system whereas the PLL assumes a wideband system. 

The local oscillator is initially assumed to be sta¬ 
bilized in frequency and amplitude, and later a frequency 
shift is introduced. The filter manages to acquire and 
track the phase in a high carrier-to-noise ratio (CNR) en¬ 
vironment although it does not estimate the frequency shift 
well. At 30 db CNR, the phase is estimated within about 
0.02 radians, whereas the frequency shift estimation error 
is on the order of 500 kHz for a frequency shift of 2MHz. 
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OPTICAL PHASE ESTIMATION 


FROM INTEGRATED SAMPLES OF THE 
HETERODYNED WAVEFRONT 

I. Introduction 

The phase of an optical wavefront cannot be measured 
directly by a detector since this would require an extra¬ 
ordinarily fast detector and associated processing elec¬ 
tronics. All presently used optical detectors respond only 
to the incident total optical power, not the instantaneous 
value of the E-field. Nevertheless, optical phase can be 
measured indirectly via several methods, most involving 
some sort of heterodyning system which moves the carrier co¬ 
herently to a lower frequency where either the spatial or 
temporal characteristics of the resulting interference pat¬ 
tern can be examined (Ref.4:3301). Among other things, op¬ 
tical phase can be used to determine refractive index vari¬ 
ations, examine transparent objects, perform nondestructive 
testing of surfaces, and reconstruction of an incident light 
wavefront <for instance, to correct aberrations in an active¬ 
ly compensated optical system, whether introduced by the op¬ 
tics or the atmosphere). 

As is demonstrated in the theoretical development later 
in this paper, the output of a detector in an optical 




receiver is at the heterodyne frequency, typically on the 
order of tens of megahertz or more. In order to pass this 
waveform, the detector must have a bandwidth of at least 
that much, and usually more if doppler shifts may be en¬ 
countered. Such wideband detectors are expensive and not 
readily available in large array structures with uniform 
response characteristics. On the other hand, large uniform 
arrays of charge coupled devices (CCD’s) are readily avail¬ 
able. These devices do not have sufficient bandwidth to 
produce an output which is directly proportional to the het¬ 
erodyned signal and thus produce a vector of outputs that 
are samples of the time integral of the heterodyned wave- 
front. That is, CCD's are photon counting detectors. 

Interestingly, Wyant has shown (Ref.20:2624) that, in 
the absence of noise, it is possible to determine the phase 
of a wavefront from these integrated samples. The purpose 
of this paper is to extend and generalize Wyant's analysis 
and to analyze how well a system performs in the presence of 
corruption by white, gaussian noise. 

In this paper, the phase of an optical field is defined 

as the argument of the received optical field phasor, 
j(27Tf o t o +0) 

e (in complex notation), minus the argument of 

j2ir(f o-f IF )t 0 

the local oscillator phasor, e , minus 

2TTfj F to . The time t 0 is the time of the measurement, 
i.e., the time elapsed since the chosen origin time, t=0 , 
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when phase was first defined. Clearly, without some sort 
of feedback between the received and local oscillator fields 
which defines a unique t=0 , (as is done in a PLL), the 

measured value of the phase will be a function of what t=0 
is chosen. However, changes in the phase, whether temporal 
or spatial, will be invariant with respect to the choice of 
the origin of the time axis. 

Overview 

This paper begins with a general analysis of how the 
heterodyne or intermediate frequency (IF) signal comes about 
and a short description of why the white, gaussian noise 
(WGN) assumption is valid for the cases of interest. Next, 
the preprocessing forced on the system by the CCD detectors 
is examined and a Cramer-Rao lower bound for the error com¬ 
mitted by any processor utilizing these samples to estimate 
phase is derived. Estimators are derived for various as¬ 
sumptions concerning the time correlation of the phase pro¬ 
cess using both maximum likelihood (ML) and state-variable 
approaches. Finally, the various processors are simulated 
and the results are presented and compared to the perform¬ 
ance of a first order linearized phase-locked loop (PLL). 
This is a standard benchmark for comparisons in phase esti¬ 
mation problems. Throughout the paper, only one detector 
in the array is examined and no spatial information is in¬ 
corporated. 
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II. Receiver and Preprocessor Theoretical Background 


In this chapter, the optical heterodyne receiver is 
very briefly described and a justification for the WGN as¬ 
sumption is presented. Next, the preprocessor structure is 
examined and a Cramer-Rao (CR) lower bound is derived for 
any estimator using such processed samples. Some limiting 
assumptions are imposed on the bandwidths of the phase pro¬ 
cess and the amplitude variations of the IF signal. Finally, 
the Cramer-Rao bound is compared to that of a processor 
utilizing the signal without integration or sampling. 

Optical Receiver Structure 

An optical heterodyne reciever is shown schematically 
in Figure 1. The received field and the local oscillator 
field are both given in complex notation with an implied 
time dependence of exp( - j 2tt f „ t) , where f 0 is the mean 

optical frequency of the received field. Each element of 
the detector array responds to the total incident optical 
power. The output of the rn,n-th detector is: 

i mn (t) = hf7/l U R (x m' y n' t) + U L (x m' y n’ t) 

A , 
d 

exp[-j2irf IF t] | 2 dx m dy n (1) 
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Figure 1. Receiver Schematic Diagram 

h77 / [l U R <x m' s 'n' t) l ! + l U L (x m > ! 'n’ t) l J+2Re{U R< x m -y„' t ) 

A d 

* U L (x m’ y 9 ' t)exp [ J2,Tf IF t ] } ] dx m dy n (2) 

h h 1 { l U R (x m' y n> t) l 2 + l U L ( V y n’ t) l 2+2 l U R (x m’ y n' t) l 
A d 

* * cos ( 2lTf TT? t+ 9( x m »y«>t)))dx m dy n (3) 


IF'' - ”' m >J, n' 


Where n is the quantum efficiency of the detector, h is 
Planck's constant, U D and U T are the received and local 
oscillator (LO) fields, respectively, measured in the detect 
tor plane, u R( x m ,y n ,t)=|U R (x m ,y n ,t)|exp[-je(x m ,y n> t)], and 
the integration is performed over A^, the area of the m,n-th 
detector. The local oscillator is assumed to be stabilized 
in frequency and assumed to produce a uniform plane wave so 
that U T (x .y t)=U . Assuming that e(x,y,t) does not vary 

over the area of an individual detector, factoring out 












/l d |U R (X m' s 'n' t)|2 + |U Ll idx m ds 'n and deflnln B ; 

‘mn ' mf, (|U R (x m' ! 'n' t) l ! + i U Ll !)dx m ds 'n 
A d 

' 2 l¥V y n- t )M u L ldx m d y n 

= { l U R ( W t) l' + f U L^ dx m dy n 
A d 


reduces equation (3) to: 

W 1 * " 1 m n[ 1+Y ° os(2,f IF t+0 mn (t) )] < 6 > 

Call y the modulation index and i the amplitude. It 

mn 

is easily seen that • Since only one detector is be¬ 

ing examined henceforth, the notation will be simplified by 

replacing i with i and e (t) with e (t) . Equa- 

mn s u mn 

tion (6) is the signal model which acts as the input to the 
preprocessor structure (integrate-and-dump filters) imposed 
by the CCD's. 

There is an inferent uncertainty in the exact arrival 
and conversion times of incident photons at any particular 
detector imposed by the statistical nature of light. These 
uncertainties are known to result in a noise process associ¬ 
ated with the signal that is Poisson distxibuted (Ref.5:53). 
However, if the arrival times of the photons are sufficiently 
close together, then for any finite bandwidth detector, the 
central limit theorem forces the resulting filtered statis¬ 
tics to be gaussian. (Ref.3:124). The smaller the bandwidth 
for a given arrival rate, the more closely the resulting pro¬ 
cess approximates a gaussian distribution. For typical local 
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oscillator power levels, the average photon arrival rate is 
so high that even relatively wideband detectors will exhibit 
gaussian statistics at their outputs (Ref.12). For the low 
bandwidth (compared to photoconductive or photovoltaic de¬ 
vices) CCD devices considered in this paper, the gaussian 
assumption is extremely good for local oscillator noise lim¬ 
ited cases. 

Preprocessor Structure 

The output of a heterodyne detector with noise corrup¬ 
tion is: 

r(t) = h(t )+n(t ) = i s [l+Y cos(2TTf If t + 0 (t))] +n( t) (7) 

where h(t) is the phase modulated carrier, 0(t) is the 
phase process to be estimated, and n(t) is a white, gaus¬ 
sian, zero mean noise process which accounts for the detec¬ 
tor noise. Since r(t) is the output of a total power de¬ 
tector, it clearly can never be negative. But, if n(t) is 
truely gaussian, then there is a finite probability that it 
will take on values such that r(t) does indeed become neg¬ 
ative. The resolution of this apparent contradiction lies in 
the fact that n(t) is never precisely gaussian. The noise 
n(t) arises as a result of the fact that the photon arrival 
times are not known exactly, even if the field is determinis¬ 
tic and known. For a field with intensity (power at the de¬ 
tector plane) I(t), the detector output will be 
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r(t) = al(t)+n(t) , where a is a constant of proportion¬ 
ality. It can be shown that the signal r(t) is a Poisson 
process with rate parameter al(t) (Ref.5:53). The mean 
of the process is given by otl(t) and so is the variance. 
The central limit theorem requires that, with a finite 
bandwidth detector, as the number of photoelectron conver¬ 
sions per second approaches infinity, the variable (r(t)- 
aI(t))//al(t) will become gaussian with zero mean and unit 
variance (Ref.3:124). That is, r(t) becomes gaussian with 
mean al(t) and variance al(t) , hence, n(t) becomes 
gaussian with mean zero and variance a I(t). However, the 
photon count per second approaches infinity only if I(t) 
becomes infinite. This would force the mean and variance 
of r(t) to infinity as well. For a large (but finite) 
photon arrival rate, the process will becomes very close to 
gaussian around the mean, but cannot be accurately modeled 
as gaussian far into the "tails". This is so since a Pois¬ 
son process is defined only for (discrete) positive argu¬ 
ments, the process r(t) , a Poisson process, is defined 
only for values of its pdf greater than zero. That is, 
fr(t)( r (t)) = 0 for r(t) < 0 . Therefore, the probability 

of n(t) taking on values less than -al(t) is zero, and 
n(t) is not truely gaussian. Since al(t) is very large, 
the truncation of the tails occurs at a very large number 
of standard deviations from the mean for even moderate in¬ 
tensities. Hense, the assumption that n(t) is zero mean 
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and gaussian is very good for large fields incident on the 
detector (many photo electron conversions per second com¬ 
pared to the bandwidth of the detector). See Davenport and 
Root (Ref.3:124) and Karp and Gagliardi (Ref.5:130). 

The effect of using a CCD detector is to introduce an 
integrate-and-dump filtering (IDF) on the signal before 
reaching any of the receiver electronics at the IF stage. 

The integration period of the CCD is user controlable as 
long as it is more than a certain minimum value set by the 
device characteristics. A typical minimum value might be 
O.lpsec (Ref.2). Designate the number of integration per¬ 
iods per period of the IF by P^ . The output of the IDF 
is a sequence of random variables with the form: 

"^ + ^ (i)+27rmi 

r i = / r(t)dt 

- gp" + p-( i-l)+2Trm( i-1) 
t t 

p-( i-i)+27rmi 

= J h(t) + n(t) dt (8) 

J-(i=|)+2Tim(i-l) 

where T=l/fj F . The factor of -T/2P t occurs in the in¬ 
tegration limits in order to take advantage of trigonometric 
simplifications that result. The factors of 2iTm are 
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included in the limits to indicate that it is not necessary 
to integrate over just one IF period. Since the IF per¬ 
iod, T , is typically on the order of 20 ns or less, and 
the CCD connot integrate for less than a few tenths of a 
usee, the CCD will obviously integrate for more than T/P 
seconds per sample. Since h(t) is 2tt- periodic, adding 
a factor of 2iTm in the limits does not change the samples 
except that the variances of the outputs of the CCD's are 
reduced by the additional integration time and the constant 

bias (due to integrating i ) is increased deterministically. 

s 

Throughout this paper, m is set to zero simply for con¬ 
venience. If m^O , the restriction on the upper limit to 
the bandwidth of the process 9(t) , to be imposed shortly, 

will be more stringent. As shown in equation (8), the se¬ 
quence {r.j.} is the sum of a signal sequence {h,^} and a 
noise sequence {n^ . The output of the detector is ob¬ 

served for a total of K IF periods, the length of the ob¬ 
servation interval is KT second and the length of the ob¬ 
servation vector is KP^ elements. If m^O , the total ob¬ 
servation interval is K(mP t +l)T . The observation vector 
is composed of K subvectors each of elements: 


£ * ( r i, r 2 


r P t :r P t +l* * *’ ,r 2P t ; *"* ;r (K-l)P t + l, 


• • • » r |£p ) ( 9 ) 
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(10) 


, nT 

= (£ i .£2 . I K } 

The parameter P has been left as an arbitrary inte¬ 
ger. If there is a value of P that minimizes the mean- 
squared error of the processor, that would be the optimum 
value of P to use in a processor. To determine if such 
a value exists and what that value is, a Cramer-Rao (CR) 
lower bound is computed. 

Cramer-Rao Lower Bound 

The CR lower bound for the unbiased estimate of a par¬ 
ameter c is given by (Ref.17:66): 

e[u — £) 2 ]>e 2 = -E[^i In f(rlc )]" 1 (11) 

where £ is the estimate of the parameter, r is the total 
observation vector, f(r|s) is the probability density func¬ 
tion (pdf) or £ conditioned on knowing the true value £ , 
and E^J is the expectation operator. The conditional pdf 
f(r|o(t)) is required. The signal is a deterministic func¬ 
tion of the phase and the noise is assumed to be independent 
of the phase. Since the noise is additive, f(r|0(t)) is 
just f(n) shifted to a mean value of h 

Since n(t) is white and gaussian, and the integral 
is a linear operator, integrated samples of n(t) will form 
a gaussian sequence. It is easily derived that {n^} is a 
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white and zero mean sequence, as well. 

E [ n i] = E [ / n ( t ) dt ] = J*E[n(t)] dt = 0 (12) 


E [ n i n j] = e[ jn(t)dtj n(t')dt'j (13) 

i j 


= yyE[n(t)n(t')]dtdt' (t-t" )dtdt" (14) 


N 0 T 


6, . 


2P t ij 


(15) 


The two-sided spectral density of n(t) has been denoted 

N 0 /2 . The appearance of the Kroenecker delta, 6.. , 

^ J 

demonstrates that the sequence is white, as expected. The 
pdf on n can now be written directly: 


f/ n \ -( 2P t Vt K 

f(n) = I - I exp 

\/2iN 0 t ' 


1 2P t 

2 *N 0 T 


Similarly, the conditional pdf of r|6(t) is: 
2P. \P^K 


f(r|0) 


. ( 2P t ¥V 

\/2ttN 0 T/ 


exp 


P t K 


2 ’ N 0 T 2 ( r i h i)‘ 
i=l 


(16) 


(17) 


From equation (11) the CR bound is: 

P. 


‘6*, t y*>, 

= -E Jq t( ~ N 0 T iL (r i -h i> > 

• ■ 


(18) 
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I 

'i 




( 



(19) 


Where the fact that ®[ r i] = was ex Pl°ited (since n(t) 

is zero mean) to remove the term involving the second deri¬ 
vative of h^ . Inserting equation (6) into equation (8) 
gives: 


h i (6(t.)) 


j. TY 
s 

2tt 


cos(pi(i-i)+0(t i ))-cos(^-(i-|)+0(t i )) 


( 20 ) 


where t i =T(i-J)/P t . It was also assumed in this integral 
that 0(t) does not change over the period of integration, 
else the integration could not be performed without explicit 
knowledge of the functional form of 0(t) , if then. Fur¬ 
thermore, if the phase changes much over one measurement in¬ 
terval, then some sort of average of the phase over the in¬ 
terval will be estimated. Consequently, it will be assumed 
that the bandwidth of the process 0(t) is small enough 
that the realization being examined changes very little over 
the measurement interval, with high probability. Equivalent¬ 
ly, some statistical measure of the bandwidth of the phase 
process is assumed known a-priori and the measurement inter¬ 
val is chosen so that the phase process realization will 
change by less than some small, predetermined bound during 
the measurement interval with some high degree of probabil¬ 
ity chosen by the designer. If the m in equation (8) is 
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I 

I 

i 


! 


t 


i 





large, this can be a rather stringent condition. 


With these assumptions, the CR bound becomes: 


= §f(rV 2 

S L t . 


X = (i-1) + 0 i 

t 


( 22 ) 


Using the trigonometric identity for cos(x)-cos(y) and apply¬ 
ing the results of Appendix A to the summation results in: 


- irCr 1 ^ 2 [ 2 p t 2 * K ' sin 2 (|r)l~ l 

S L ** j 

Tn 0 2 1 it 2 

[2TK * (i s Y ) 2 J * sin 2 (TT/P t ) 


N o 

As P^*® , this equation reduces to £2= 2 gr*' (i y) 2 which 

s 

is the same as the error for the linearized first order PLL 
(Ref.18:93), or the CR bound for a system which has access 
to the non-integrated, non-sampled waveform. For P^ fin¬ 
ite, the error bound will be greater than that for the PLL 
by a factor of tt 2 /P 2 sin 2 (p—) . In all simulations, P t = 

4 simply to reduce computer time requirements. For larger 
P , an appropriate scale factor can be introduced into 

v 

the simulated errors by using equation (24). 
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III. Processor Determination : 
Maximum Likelihood Estimation 


Three specific cases of the estimation problem can be 
identified: 

(i) 0(t) is a constant, 8 0 , for all time. It may 

be construed as a single realization of a random variable 
uniformly distributed between -n and tt radians. 

(ii) 0(t) is a random process with unknown statistics 

(iii) 0(t) is a random process with known statistics. 
The first two cases are analyzed in this chapter and it is 
shown that the optimum processor is a discretized phase-lock 
ed loop. Because of the nature of the integrate-and-dump 
preprocessing imposed on the measurements, it may not be ob¬ 
vious that the optimum ML processor still has a PLL struc¬ 
ture. The third case is analyzed in the next chapter. 


Measurement Model Equations 

Before continuing, it is convenient to derive the exact 
form of h^(0) and its first two derivatives with respect 
to 0 . Using equation (6) in equation (8) gives Ik(0) , 

which can be differentiated directly. 


h^Q) 


igT jr + 7 sin(J-)cos(p^-(i-l) + 0) 

_ • 
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6h (6) i TY 2 

-TT- = - -T- sin(^)sin(|-(i-l) + 8) 


6 2 h (6) i TY 2 

= _ -f_ sin( l_ )cos( |l (i _i )+ e) 


These equations are used many times in the derivations that 
follow and will be referred to quite often. These equations 
are needed because the first derivative is always necessary 
in the derivation of the ML processor. The second derivative 
will be used when the Kalman filter is derived because the 
nonlinear measurements will be linearized in a Taylor series 
expansion and the second derivative appears in the deriva¬ 
tion . 


Maximum Likelihood Approach: Case (i) 

The ML estimate of 0 O is the value of 0 O that maxi¬ 


mizes the conditional density f(r|0 o ) . Since f(rj0 o ) 


is gaussian, the ML estimate is the value of 0 


o . o o > 


that solves: 


■jg— In f(rj0o) - 0 


Taking the logarithm of equation (17) and differentiating 


gives: 


V' ° h i 

^ < r i- h i>607 = ° 
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Inserting equations (26) and (27) results, after some trig¬ 
onometric identities are applied: 


0 = ^r. sin(p-(i-l) + 0 o ) - ^ 


i T 

p^- sin(^-)sin(|^(i-l)+( 
t t t 



Yi T 
s 

2tt 


sin(p-) 

t 


sin(il(i-l)+2G 0 ) 

t 


(30) 


The trigonometric terms in the last two summations are 
shown to sum to zero in Appendix A. The estimator reduces 
to: 

KP t 

£ r. sin(|l(i-l)+0) = 0 (31) 

i=l z 

This processor is a discrete analog to the familiar 
PLL and can be implemented in closed loop form for any val¬ 
ues of K and P^. . A block diagram is shown in Figure 
2. The dotted box encloses a discrete IPF(DIDF). 

Maximum Likelihood Approach: Case (ii) 

Suppose that 0(t) is not constant over the entire ob¬ 
servation interval. In light of Appendix A, let the phase 
be approximately constant over one period, T , but vary from 
period to period in an unknown manner. Equation (31) be¬ 
comes : 
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If equation (32) is to be true for any value of K , then 
the inner sum must be zero for every m . The estimator 


is then: 


E r i sin<§;<i-l> + V - 0 

i=l x 


The estimator produces a sequence of estimates, {§ m } • As 
might be expected, 0 m is estimated independently of any 

/s 

past value of 0 or 0 and only data from the current 
mm 

measurement interval ^mT, (m+l)Tj is used to perform the 

T 

estimation. That is, r=(r ,r ,...,Tp ) . This is direct¬ 
ly a result of the fact that {6.} is (or is assumed to be) 


an uncorrelated sequence: 


0 kl ' 
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IV. Processor Determination: State Variable 


Analysis and Kalman Filtering 


If some information concerning the time correlation of 
the process 0(t) is already known, it should be possible 
to obtain a better estimate (lower variance estimate) of 
the process by incorporating the a-priori information into 
the processor. In this chapter, state variable models and 
nonlinear Kalman filtering are applied to the problem of 
incorporating a-priori information into the processing al¬ 
gorithm. A few words are said about state variable equa¬ 
tions in general and then the arguments are specialized 
for the measurement model of interest. 

State Variable Equations 

The phase process is assumed to be the output of a lin¬ 
ear, time invariant system driven by white gaussian noise. 
This system is described by n state variables in a vector 
differential equation (Ref.13:556). 


x(t) = F x(t) + Gu(t) 

(34) 

9(t) = c T x(t) 

(35) 

c T = (0,0,...,0,1) 

(36) 


where: 
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x(u) n-dimensional state vector 
u(t) m-dimensional noise vector 
F nxn-dimensional system matrix 

G nxm-dimensional noise system matrix 
The solution to equation (34) is (Ref.13:556): 


x(t) = 4>( t-t o )x( t a ) + / f|)(t-t^)Gu(t')dt' 


Where 4> (t—1 0 ) is called the system transition matrix and 
solves equation (38) with initial condition (39)(Ref.13:557) 


<Kt-t 0 ) = F(J)(t-t o) 


<K 0) = I 


Since the data to be processed is discrete, a discrete for¬ 
mulation of equation (37) is: 


x(kT) = x k = [<KT)] k xo + 2 [<K T )] k u 1 (40) 


This equation can be put into a recursive form with only mod¬ 
erate algebraic manipulation (Ref,13:557): 


2k = *2 k -l + 


For the time invariant case, the solution to equation (38) is 
exp(F( t-t o) )=<J)( t-t o) which becomes, for discrete measure¬ 
ments (Ref.13:558): 
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<f> = <ji(T) = exp^FTj 


(42) 


It will now be assumed that the phase does not change 
during an IF period (or changes very little) so that the 
phase process is described by a sequence of random variables 
{0^} or, alternatively, {x^} , where the variables are sep¬ 

arated by T seconds. This will allow taking advantage of 
the trigonometric simplifications proven in Appendix A. 

-k = - (6 k ) + -k = -^£ T ik ) + -k (43) 

Where n k is a white, gaussian vector with statistics: 



The dimensionality of r, h, and n is . The matrix 

I is a P^xP^ identity matrix. The sequence {u^} is 
white and gaussian with covariance matrix Q . The sequences 
{n^} and {u k } are independent. These assumptions, plus the 
initial conditions are: 


E 


■fej' 2 

E [n k 3ll] - « 5 kl 

(45) 

r Ti 

r Ti 


Uk-lJ “ 0 

E [—k—1J = 

(46) 

E[x : ] = x, 

E [(x-x 0 )(X-Xo) T ] = Po 

(47) 
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Given the received sequence {r k l . an "optimum" es- 

A 

timate of the state {x^} is desired from which an optimum 
estimate of the phase can be derived, or any 

other element of the state vector. Unfortunately, no gen¬ 
eral solution to this problem exists for arbitrary nonlin¬ 
ear functions h(•) in either the continuous or discrete 
case. Exact solutions do exist for linear h(•) and ap¬ 
proximate solutions do exist for a nonlinear li(• ) which 
has been linearized by a Taylor series expansion about the 

A 

best estimate of the state, x^ . Sage solves the prob¬ 
lem using a discrete invariant embedding procedure (Ref.14: 
450) to give the equations for a maximum a-posteriori (MAP) 
estimate: 

I k = h ( 2 £ k )+n k (48) 


—k+l " + Vi' 


-T 6 U T 


Sx k 


h ()H"^ k+1 -ii()] (49) 


k+l 


’ -T a -T 

6 

<J> -<p 

6^T 


~ k . 


-4- n T u£ k > R 1 
% 


r k+1 -h(<J>x k ) 


r lp k 

(50)' 


p k = ^ P k + g Q qT ^ T (5i> 

A T A A 

When h(4 ) Xj c ) == ll(c (})X^)=h( 0^) , this can reduce to (see Appen- 

dix B): 

-k = ^ (6 k ) + Hk (52) 
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where: 


CNR 


( V>7 N. 

2 / 2T 


(62) 


Z k+l (l) 


i Y 
s 


sin(p-)cos(|i(i-l)+0 k+1 )+n^ +1 (i) 


(63) 


2 = No 

n" 2P t T 


(64) 


The carrier-to-noise ratio (CNR) thus defined is the ratio 
of the signal power to the noise power in the total band¬ 
width used in producing the estimate (1/T) . This is the 

same CNR used by Viterbi (Ref.18:93) and thus allows direct 
comparisons with his results for the PLL, In addition, the 
CNR gives an intuitive feel for how "good" the measurements 
are, or how badly they are corrupted by measurement noise. 


Processor Ambiguity 

Appendix B shows that the processor equations can be 
further simplified to: 


—k+1 



P k+1 -[ B ’ sln(0 k+l“®k ) + n ks] 


(65) 


p k :i ' K* T+ «']" + ££ T [6' 


c °s<9 k+ i-e'Hn kc 


( 66 ) 
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where: 


P 2 

6 = ^-r* sin 2 (p—) *CNR 


(67) 


var 



= var 



= 6 



( 68 ) 


The sequences {n kg } and (n k } are also zero mean. For 

any plant model which is one dimensional, —k~ x k = ^k anc * 

T 

cc =1 . For this case, plugging equation (66) into (65), 

assuming the noise is negligible, and dividing top and bot- 
tom by cos(0 k+1 ~0^) : 


9 k + l = 6 k " 


tan(0 k+ i-9 k ) 


8"^P k ct) T +Q")"sec(e k+1 -e') + l 


(69) 


For large CNR, the estimate should be quite good and, with- 
in modulo-2TT, 0 k+l -0 k~^ ' Since sec(0)~l , equation (69) 

reduces to, approximately: 


6 k + l " 0 k - 


- 1 


1+8 (*P k f<) 


-1 




(70) 


However, the function tan(») exhibits a modulo-n ambiguity, 
so that tan(e k+1 ~( § k +pTT) )=tan( e k+1 -e k ) , where p is any 
integer. Thus, the estimator will be unable to determine 
the phase better than modulo-n 

The above arguments were for a one-dimensional state 
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vector with high CNR. Simulations have confirmed this re¬ 
sult for smaller CNR's as well. 


Constant Phase 

The case of a constant phase plant is the simplest 
case to model. The state vector has but one component and 
the matrices collapse to scalars. The state variable equa¬ 
tions become: 


9 


k+l 


(71) 


4> = 1 G = 0 R = Q = arbitrary (72) 

The estimator or Kalman filter equations are, from (65) and 
( 66 ): 


J k+1 


= 9, - P 


k+l 


[e.sin(0-0 k )+n ks ] 


(73) 



Figure (3b) is a block diagram of this processor. Figure 
(3a) shows the processor in the form given by equations (57) 
and (58). 

Random Phase Walk 

The random walk is modelled by the state variable equa¬ 
tion : 

9 k+l = °k + Gu k < 75 > 
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Figure 3a. Constant Phase Processor Structure (Realization 1), 
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]=Q=1 , 
of the noise term is 
$=1 , and R=N 0 T/2P t 
filter equations are: 


Let 


var 


K 


for simplicity, so that the variance 
reflected in the value of G . Also, 
as in the constant phase case. The 


k+1 


= 0, - P 


k+ 


it 


S-sin(6 k+1 -e k >n ks j 


(76) 



(P k+ G 2 ) 1 + Bcos(0 k+1 -0 k ) + n kc 


(77) 


Figure (4b) shows this processor while Figure (4a) shows 
the implementation according to equations (57) and (58). 
Both forms of the filter are shown because while the form 
of equations (65) and (66) gives a good intuitive feel for 
how the filters process the estimation error, the form of 
equations (57) and (58) is the actual form a hardware pro¬ 
cessor would take. This is true because actual phase is 
quite clearly unavailable to the processor. 


Demodulation of FM Signals 

Suppose that the information of interest in a state 

T 

vector is given by m(t) = c x(t) • Suppose also that the 
phase of the carrier is given not by m(t) , but by 0(t)= 

t . t 

X/m(i)dT + 0o so that 0(t)=m(t)=Ac x(t) w'here A is called 
the frequency modulation index. A new state vector can be 
constructed from the old state vector and 0(t) . The new 

state variable equation is: 
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X (t) 


x(t) 


' F 0 ' 


"x(t) 


"g 0 ' 

u(t) 


= 




+ 



_0(t)_ 


1 - 

10 

h- 

*4 

CD 

1 _ 


0(t)_ 


-° G d 

u Q (t)_ 


x'(t) = x'(t) + G' u'(t) (79) 

Where F. has been added to the F' matrix to account for 
possible variations in 0(t) not due to the FMing of the 
carrier by m(t). In the discrete formulation: 

2£k+i = + 4>' G 'Hk (80) 


<T 




(81) 


The plant equation, and hence the filter equation, is at 
least two dimensional, but this is no theoretical compli¬ 
cation. 


A practical case is when the frequency is a constant 
over the observation interval. This could come up when 
estimating phase in the presence of an unknown doppler 
shift. Model the system with two state variables and sup¬ 
pose the phase, in the absence of the doppler shift, would 
be executing a random walk. Then: 


r 

o 

o 


x(t ) 


i 

o 

o 

_1 


'u(t) 




+ 




I 

o 

«-< 

_ 1 


1 

V/ 

O 

_ 1 


1 

o 

a 

CD 


- u e (t )- 


( 82 ) 
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= (4,^P k > T +G'Q-G' T ) ‘+ c c T Bcos(0 k+1 -0')+n kc 


where: 

E[u'u'] = Q"6 kl 8£ = cV£ k (87) 

It is now possible to estimate both phase and frequency. 

rn 

The doppler shift is given by x k+1 =d x k+1 and the phase is 
given by \ +1 = £ T E^+l » where d T =(l,0) 

Phase Estimation with Frequency Uncertainty 

Imagine a case where it is desired to estimate the phase 
of an optical field, but the local oscillator is known to be 
unstable. Alternatively, the LO may be stable and the re¬ 
ceived field frequency fixed, but the integration limits may 
not be controlled precisely. The uncertainty in the limit 

times, =—(i-£) , can be modelled as an uncertainty in f T = 

Ir 

1 
T 
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The frequency uncertainty can be incorporated into the model 
of equation (80) by simply allowing G=^0 . Unfortunately, 

the state-variable models assume a driving force which is 
white, and the frequency uncertainty of a local oscillator 
is not even nearly white for any case of interest. The 
state-variable models can accomodate non-white driving func¬ 
tions by driving a "coloring filter" with white noise and 
allowing the output of the coloring filter to drive the sys¬ 
tem of interest. The state variable equations for the col¬ 
oring filter are then used to augment the original state 
variable equation (Ref.7:180). The coloring filter is de¬ 
scribed by: 

x f (t) = F f x f (t) + G f u f (t) (88) 

u(t) = H f x f (t) (89) 

Where the notation in equation (89) is meant to indicate 
that only the u(t) portion of u'(t) in equation (79) 
is augmented while u Q (t) (and, for that matter, any por- 
tions of u(t) that are already white) is not an output 
of the system x^(t) • The augmented system is described 

by: 
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( 90 ) 


( 91 ) 


The 


x(t) = F • x (t) + G 
—a a —a 

discretized equation is 


• u (t) 
a —a 


given by: 


x (k+1) = 
—<x 


i x (k)+if> G u (k) 


(92 ) 


*a = exp[F a T] (92) 

Where, for ease of reading, the subscript k has been re¬ 
placed by the argument k . The filter equations are sim¬ 
ply equations (65) and (66), or equations (57) and (58), 
with all quantities replaced by their associated augmented 
counterparts. Note carefully that u (t) is now a white 

3 . 

process driving the augmented system. This is as far as it 
is possible to go without assuming certain correlations 
for x f (t) 

It is cautioned that it was assumed that G(t) could 
be modelled as the output of a system driven by white noise. 
If this is not the case, then equation (89) is modified so 
that u (t)=H^x f (t) . Also, it was assumed that u(t) is 

the output of a system all of whose inputs were non-white. 

If this is not the case, then some elements of u(t) are 
not included in (89), are not augmented, and enter equation 
(90) in their original form. A further caution is that 
since the white noise driving function, u^(t) , is gaus- 

sian, the resulting density for u(t) , the frequency 


36 





uncertainty, is also gaussian for any linear plant model. 
This is because a gaussian process remains gaussian after 
being transformed by a linear system. So, if a white, gaus¬ 
sian process u^(t) drives a linear, time invariant system 
to produce an output u(t) and u(t) drives a laser mod¬ 
eled as a voltage controlled oscillator centered at f 0 , 
the pdf of the laser output frequency is the same as the pdf 
of u(t) , i.e., gaussian. A nonlinear plant would be re¬ 
quired in order to accurately model a homogeneously broad¬ 
ened laser line, which is Lorentzian, since only through a 
nonlinear transformation can a gaussian process u f (t) be 
transformed into a nongaussian process u(t) . A Lorentzian 
random variable is also known as a Cauchy random variable. 
The pdf of a Cauchy distributed random variable is given by 
(Ref.19:48): 


f 


x 



_ a 

a z +(x-m ) 1 
v x 


(94) 


(Technically, of course, none of the moments of this density 
exist. That is, /x n f x (x)dx does not converge over infin¬ 
ite limits for n>0 ). Of course, an approximation would 
still be made later, when the estimator equations are lin¬ 
earized in a Taylor series expansion and this will result 
in a different modeling error. 
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V. Simulation of Kalman F jIter Processor 

In this chapter the signal models of chapter IV are 
simulated and the results of the simulations are discussed. 
Comparisons are made with CR bounds and a linearized first 
order PLL. Some interesting behavior on the part of the 
processor is noted and explained. 

Simulation Parameters 

A subroutine was written in Fortran IV which performs 
a Monte Carlo simulation of the processor of interest. The 
simulator performs an ensemble average over a user deter¬ 
mined number of realizations of the measurement noise pro¬ 
cess. The subroutine is listed in Appendix C. 

For the case where the IF is fixed and known (no dop- 
pler shift is present), the phase was modelled in one di¬ 
mension and two models were investigated. The first model 
was for the phase a constant. The second model was for 
the phase executing a random walk around the unit circle. 
These two models were chosen because they represent two ex¬ 
tremes in the degree of correlation of the phase process 
and are typically the models encountered in the literature 
(Ref.13:939, Ref.8:48). Neither the constant phase plant 
model nor the random walk phase plant model may be adequate 
for any particular case, but since they bracket the process 


38 








correlations likely to be encountered, they can serve as 
upper and lower bounds, respectively, on how well the fil¬ 
ter will perform in practice. 

Also simulated was the case where the phase is execut¬ 
ing a random walk and there is a constant, unknown, doppler 
shift. The errors in jointly estimating both phase and fre¬ 
quency are investigated. 

In all simulations 100 sequential time estimates of 
the phase process are computed and each estimate is en¬ 
semble averaged over 100 realizations of the noise process. 
One hundred realizations of the noise is not very many, but 
it is common practice in Monte Carlo simulations to use a 
number of realizations that is relatively small compared to 
what the central limit theorem might deem prudent for a 
reasonable confidence in the statistics (Ref.16:939, 15:382, 
10.67). Compared to the references cited, 100 is a rather 
large number of realizations. The reason that so few real¬ 
izations are used is the quantity of computer time consumed. 
One run of the subroutine with 100 realizations and 100 re¬ 
cursions will consume about 10 to 15 seconds of CP time. 

For a 0.9 probability that the phase error estimate is with¬ 
in 0.1 radians of the actual phase error, the central limit 
theorem would require 809 realizations. The weak law of 
large numbers would require 3290 realizations (Ref.10:67). 
The Chernoff bound requires 1515 realizations (Ref.19:97). 
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Constant Phase Processor 

For the constant phase processor, the phase was mod¬ 
elled as a random variable uniformly distributed between 
- 7 t and 7 t radians. Initial values were then x 0 = 0’o=0 , 

Po = tt 2 /3. Po is the variance of the a-priori estimate 0 O . 
The state vector is one-dimensional. A value of P^=4 
integrals per IF period T was used. 

Figure 5 shows the simulation results for an actual 
angle of 0 =0 and three different CNR’s of +30db, +10db, 
and -lOdb. For CNR=30db, the processor converges almost 
immediately to the correct value. For CNR=-10db , the 
processor fails to converge at all after 100 recursions. 
Because the derivation of the Kalman filter performed by 
Sage discarded all terms in a Taylor series expansion of 

A 


best estimate of x^+i before measurement 


x^ and P^ beyond the first order, the estimator does not 
perform well in a high noise environment. In a high noise 
environment, the assumption of approximate linearity of the 
processor around a best a-priori estimate trajectory (the 

is incorp¬ 
orated, 2£k+i = ^2Sk ) is violated and consequently the filter 

no longer adequately extracts information from the measure¬ 
ments. As a result, the processor fails to converge. In 
the low CNR environment, the processor variance does not 
settle to a steady state value as it does in the high CNR 
environment, but rather it oscillates about, never reaching 
a steady state. The processor simply tracks the noise. 
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For CNR=10db, the processor does converge after about 25 
recursions, but it converges to the wrong value, about 
^100 = ~ 0,01 rac ^ ans (-0.57°). The processor is biased. 

Many texts either assume that nonlinear Kalman filters are 
unbiased or state without proof, justification or condi¬ 
tioning that the bias is zero (Ref.14:433). As discussed 
in Maybeck (Ref.8), these claims are false. Maybeck dis¬ 
cusses methods for reducing the bias of the processor by 
using second order approximation filters (retaining Taylor 
series terms up to order 2), or by using first order filters 
(like this one) but adding bias correction terms taken from 
the second order filter equations. These methods result in 
somewhat more complicated filter algorithms and were not 
used in this research, however, processor biases will be 
pointed out as they are encountered. 

To investigate the bias problem further, Figure 6 plots 

/N 

the estimation error after 100 recursions, ®ioO~® ' 
should be very close to the actual bias, for three CNR's. 
First, notice that for CNR=-10 db, the bias varies linearly 
(approximately) with the actual angle 0 . This is because 

the processor is tracking the noise which is zero mean and 
tends to guess 0=0 , the a-priori estimate, on the average 

Also notice that even for a very high CNR, 30db, the process 

A 

or guesses 0=0 , on the average, whenever 0=>±n/2 . This 

phenomenon is a direct result of the modulo-ir ambiguity of 
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the processor. If the actual angle is +tf / 2 , the process¬ 

or may determine an estimate that, due to noise corruption, 
is 0±e , where £>0 is the error. But |+e(or -^-e) is 
outside the ambiguity range of the processor and so it 

X "ft 

guesses -~+e (or^- -£ ) • That is, the processor "wraps 

around" the semicircle and the average over many realizations 

is zero. This explains why, even though the estimate seems 

poor from the examination of a sample mean, the variance, P , 

K 

still gets small as k gets large. It also explains why the 
bias rises, for CNR=10db, as 0 -*— . With the higher noise 

level, as the angle gets closer to - ^ , the processor will 

"wrap around" more often and the bias will rise. Similarly, 
for a fixed angle, the processor will "wrap around" more often 
for a smaller CNR, so the bias increases as CNR decreases. 

_ i _ 2 

Figure 7 is a plot of the inverse variance ( p ioO =a ^ ) 
versus the CNR. Also plotted is the theoretical CR bound x0 ° 
and the variance for the linearized PLL. The processor sig¬ 
nificantly underperforms the PLL, although this gap would 
narrow if P^ was increased. The processor almost reaches 
its CR bound with equality, however. After 100 recursions, 
the estimator is "almost efficient". 

Random Phase Walk 

The next model simulated is for the phase executing a 
random v/alk around the unit circle. For this case, the 
variance of the walk, or diffusion, parameter was set at 
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igure 7. Constant Thase Processor Variances 













Q=0.005 so that the size of the "steps" in the walk would 
be small enough that the phase can be assumed approximate¬ 
ly constant over each IF period. For this value of Q , 
the magnitude of is less than 0.1 radians (5.7°) 

with probability 0.84. Only one realization of the phase 
process was examined and the processor's performance was, 
again, averaged over 100 noise process realizations. Ideal¬ 
ly, the processor's performance should be averaged over 
many realizations of the phase process as well as averag¬ 
ing over many realizations of the noise process. This 
averaging was not performed because of the amount of compu¬ 
ter time that it would have consumed. 

Figure 8 shows the error performance of the processor 
plotted for 100 recursions at 3 CNIt's. The quantity plot- 

A 

ted is the estimator error, 9^-9^ • The reduction in per¬ 

formance at progressively lower CNR's is self evident. Note, 
however, that the processor, at lOdb CNR, tends to made more 
negative errors than positive errors. This is another in¬ 
dication of the processor bias. The bias is even more evi¬ 
dent for the CNR=-10db case. For this particular realiza¬ 
tion of the walk, the phase managed to walk as far as +1.4 
radians, which is quite close to +— . For a walk realiza¬ 
tion which does not walk as far from 0=0 , the error, on 

the average, would be expected to be less. Simulations con¬ 
firm that this is true, although the processor still makes 
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significantly more errors for the lower CNR. 

Figure 9 plots the inverse error variance versus the 
CNR along with the CR bound and the theoretical results 
for the linearized PLL. The scale on the vertical axis 
has been changed relative to Figure 7. This change re¬ 
flects the increased bandwidth of the process and the cor¬ 
responding increase in the bandwidth of the PLL loop fil¬ 
ter. The loop filter time constant has gone from KT, where 
K=k =100, to T . The fact that the processor seems to 

outperform its CR bound can be explained on the basis of 
the fact that the computation of the CR bound in equation 
(24) set K=1 . That is, the CR bound is for a processor 
which does not take into account the a-priori information 
on the correlation of the phase process. The processor 
simulated here does use this a-priori information and it 
gives a lower variance estimate as a result. 

Since the processor is biased, and the CR bound was 
derived for an unbiased estimator, it might be argued that 
the processor outperforms its CR bound because the bound 
was derived for the unbiased case. The CR bound for a bias 
ed estimator is given by multiplying the right hand side of 
equation (11), (the bound for an unbiased estimator) by 
(l-db(0)/d0 ) ? , where b(0) is the processor bias as a 

function of 0 (Ref.17:146). From Figure 6, for moderate 

to high CNR's, db(0)/dG=O for all angles not close to ± 
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At low CNR's, the Kalman filter also outperforms the PLL. 
This is also due to the fact that the PLL does not incorp¬ 
orate a-priori information about the plant dynamics, where¬ 
as the Kalman filter does incorporate this information. 


Phase Estimation with Doppler Shift 

It was mentioned earlier that phase and frequency can 
be estimated jointly. The state vector becomes two dimen¬ 
sional with components Xj=doppler shift and x 2 =phase. The 
phase process was modelled as a random walk with diffusion 
parameter Q=0.005 . A static model of the doppler shift 

was assumed with a shift of +2MHz. The IF was chosen at 
l/T=50MHz. If the filter is in a 1.06y Nd:YAG laser receiv¬ 
er, this corresponds to a relative target velocity of about 
5 mph (2.1 mtrs/sec), a rather small velocity and doppler 
shift. The following model parameters were used for the 
simulation: 



r 

o 

o 

i_ 


l 

o 

tH 


~3.6x10 1 3 1.09x10 6 " 

Q = 


<t> = 


Po = 



_0 0.005. 


_T 1_ 


_1.09xl0 6 7T 2 / 3 _ 




( 95 ) 
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One hundred recursive estimates were simulated and averaged 
over 20 realizations of the measurement noise. 

The P 0 matrix was chosen iteratively in an effort to 
get the filter to converge without tracking the noise (Ref. 
7:338). Initially, Poii was set to 36 THz 2 by assuming 
that the frequency error or doppler is normally distributed 
with standard deviation 6 MHz. The phase variance, P 022 , 
was again chosen at tt 2 / 3 . The initial trial of P 0 i 2 = Po 2 1 
was /p777p 77T=1.09x10 7 but this value was too large and 
caused the estimator to track the measurement noise. Ex¬ 
amination of the estimator equation (57) indicates that the 
element P^ 12 influences the update of xi , so P 012 
was reduced until the filter converged more quickly to the 
actual variance of its estimates (or something close to it 
in this case). 

At 30 and lOdb CNR, the filter makes the following av¬ 
erage errors and standard deviations: 

", 56MHz (o 20 )n = . 357MHz 

CNR=30 db: x 20 -x = (96) 

-.024rad_ ( 020)22 = .102rad 


CNR=10 db: x 20 -x = 



( 020)11 = .375MHz 

(97) 

( 020)22 = .159rad 


The processor exhibits little degradation in performance 
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when the CNR drops from 30db to lOdb compared to a similar 
drop when the frequency is known exactly. However, the var¬ 
iances of the errors are significantly higher than the error 
variances for the same CNR's when the frequency is known 
exactly. No comparisons were made with other processors, 
but the phase estimate seems better than might intuitively 
be expected for the given error in the estimate of the fre¬ 
quency error (^30%). This is not that surprising since the 
actual frequency is 52MIiz, an error in the estimate of the 
offset of 500KHz gives a fractional frequency error of less 
than 1%. 
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VI. Conclusions and Recommendations 


A nonlinear Kalman.filter has been derived to estimate 
the phase of a carrier measured through an IDF. Two models 
of the phase process were examined and one case involving 
doppler shift was also examined, though not in depth. The 
problem was motivated by an optical heterodyne receiver 
structure, but the solution formulation was much more gen¬ 
eral in nature. Any system which processes samples like 
those in equation (25) could use all or part of the theo¬ 
retical development in this thesis. 

Conclusions 

This study has demonstrated that it is possible to es¬ 
timate the phase of an optical field with detectors which 
are much lower in bandwidth than the heterodyned waveform 
at the IF. The only requirement placed on the system is 
that the integration time of the detector be small enough 
that the phase and amplitude can be treated as constants 
over the integration period. 

Although the processor shows no significant improve¬ 
ment over the linearized PLL for the two phase models ex¬ 
amined, the PLL could be utilized only in conjunction with 
wideband detectors. Thus, this system allows phase measure¬ 
ments where they would otherwise be impossible. Also, for 
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cases where distributed phase measurements are desirable, 
wideband detectors are generally not available in large ar¬ 
rays, but CCD's are. 

The few simulations performed have also indicated that 
it is possible to estimate relative frequency, doppler 
or LO instabilities, in addition to relative phase. These 
simulations have also indicated that the phase can be track¬ 
ed even with a relatively poor estimate of the frequency 
available. The system is not too sensitive to the actual 
LO frequency. 

Recommendations 

The measurement equation (25) was derived on the basis 
of two simplifying assumptions: the amplitude and the phase 
are both approximately constant over the measurement inter¬ 
val. Since a CCD inherently averages over a rather long 
period of time compared to T , the first assumption would 
probably be violated in a rapidly fading environment or 
where the carrier has incidental AM of wide bandwidth com¬ 
pared to the integration time. Similarly, the second as¬ 
sumption could be investigated by allowing the phase to 
drift during the measurement interval. 

The presence of a modulo-n ambiguity in the processor 
is a drawback, but it might be possible to eliminate it. 

One possible method to accomplish this was suggested by 
Meer (Ref.10:106). That method consists of checking to see 
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if 0^ is close to ± ^ • If it is, and 0^+2 

_ A 

close to i 2 and has sign opposite that of 9^ 
new estimate of ®k +2 can be d e fi ne< l: 


is also 
, then a 


0 k + l = 9 k + l ±7T 


(98) 


The sign in front of the it depends on whether 0^ was 
positive or negative. In the absence of noise, this pro¬ 
cedure could regenerate the phase sequence with no ambigu¬ 
ity. In the presence of noise, the number of errors com- 

7T i * 1 

mitted would depend on how small the parameter ^ = 2' 0 k' 
is permitted to become and how small 6 ' = \~ I 9 k+iI mu st 

A 

also "ecome before a correction is made, provided 9 
changes sign. Clearly, more errors would also be made as 
the CNR is decreased. 

This procedure is somewhat ad-hoc. A more mathemati¬ 
cally precise method would be to take advantage of both the 
in phase and quadrature information contained in the measure¬ 
ment. The update portion of the covariance matrix, the term 


Bcos( > +n kc in e 9 uat i° n (®6), or the algebraically 

equivalent part of equation (58), is the part that incorpor¬ 
ated new measurement information into the covariance equa¬ 
tion. If l®k+l -0 k^ >Tr ’ a wra P ar ound will occur and 

the update portion of the variance equation will be negative, 

A 

since cos(x)<0 when tt < | x j < 2rr . Adding tt to if 

it falls into the fourth quadrant (or subtracting tt from 
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if it falls into the first quadrant) when the update 
term is negative will regenerate the phase sequence. If the 
noise variance becomes too large, clearly errors may be made, 
but the estimator does not perform well in very low CNR en¬ 
vironments, anyway. 

Much work remains open on the problem of estimating a 
multi-parameter state vector. In particular, the most in¬ 
teresting case from an optical communication point of view 
is the unstable LO tracking a random phase walk. The fil¬ 
ter equations were derived in (92) and (93) (implicitly, at 
least), but performance simulations were only begun u’ith a 
small, fixed unknown frequency shift. Extensions of this 
■work would have practical applications to any laboratory 
or field system since LO instabilities are always present. 

Since the main advantage of CCD detectors is their 
availability in large arrays which allow spatially distri¬ 
buted measurements, it would be logical to incorporate 
spatial information into the processor. If the phase of 
the received waveform or the background noise has a known 
or assumed correlation, then the entire array of measure¬ 
ments can be utilized to estimate the phase at each pixel. 

The incorporation of more a-priori information should re¬ 
sult in a lower variance processor. Because of the conceiv¬ 
ably large size of the array, either a piecewise correlation 
over small parts of the array would be assumed, or the array 
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could be processed in a spatially recursive fashion, 
second alternative is especially intriguing since the 
tion of the permissibility of recursive processing in 
dimensions is not quite closed (Ref.11:935). 


The 

ques- 

two 
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Appendix A 

Proof that Samples of a Sine Wave Sum to Zero 


Numerous times in the text, a great deal of mathematical 

simplification resulted when terms of the form sin(^-(i-l) 

4ir t 

+ 0(t)) and sin(p—( i~l) + 0(t ) ) summed from i = l to i=P 

were set to zero. Here this identity is established. Con¬ 
sider a slightly more general case where k is an arbitrary 
integer: 


^ sin(|— k-(i-l)+0(t)) = ^ sin(2iri-^- + 0 (t)) 

i=l 1 i=l t t 

(99) 

Euler's identity states that exp(jx)=cos(x)+jsin(x ) , so, 

if it can bo established that: 


E 

i = l 


exp 


j{2ni 


k_ 

V 


2':k 


0 (t ) } 


= 0 


( 100 ) 


then the proof is completed. The term 
cap be moved out of the summation since 
on the index, i . Define a parameter 
that a geometric series in L results, 
as: 


exp 


+ J0(t) 


it does not depend 
L=exp(,j2:rk/P^ ) so 


A geometric series 


sums 
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( 101 ) 


t p 

£ 

i=l 

P t 

Since k is an integer, L =0 unless k=0 , so the sum 

in equation (100) is indeed zero for all k^O , including 
the special cases of interest, k=l and k=2 
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Appendix B 

Simplification of Estimator Equations 


The 

to 

For 


Several vector differentiation rules must be defined. 

(5 

operator is a row vector which is always multiplied 

the right of the vector on which it operates (Ref.1:983). 
a 3-dimensional vector: 

6f (x) 


f i (x) 
f 2 (x) 

f 3 (x) 


6 x i 6 x 2 


6xi 


( 102 ) 


Similarly, the operator ^ is a column vector which al¬ 
ways operates to the left of a transposed column vector (Ref. 


1:983): 


6f T (x) 

6x 


6 

6x i 



[fl(x) f 2 (x) f 3(x)J 


(103) 


Both differentiations result in matrices with elements like 
6f i (x) 

6x . 

J 

From the chain rules for these operators (Ref.1:985), 
the partial derivatives in equations (49) and (50) can be ex¬ 
panded : 
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( 104 ) 


s , T , T, - . 5 <= T *i k > T «i> T <°£> 

h (£ <t> 2£ k ) = --- * -~- 

-k k fix,, 60,; 


T fih (ep 

5, c --- 

&&' 


(105) 


6 _T 


~T H(£ 0 * k > 
6 *k 


fih(8p T 


6o; 

k 


(106) 


Equation (49) becomes: 


-T x ^ (6 k ) „-i 

60, 


^k+i = ^k +p k+i^ —rf- R ‘ r k +i"£ ( £ ^k } 




(107) 


6 r( 0 P „.i 


-k+1 *-k + P k+1— R —k+1 - (0 p] 

60 k L J 


(108) 


which is equation (53). Equation (50) becomes: 


P = x>' 
k+ 1 k 


-T -T 

<P ~<P 

6 r*T ,S - T(0 k ) R - 1 ] 

ik+i-^^P 

- 

T , 

r> J P 



- P k 


60 k L 60 k 




- i 


(109) 


= p; 


-T 


-c-M 

L 


fp A 

6h 1 (0p 

-— R 


, l 6e; 
k L k 


Ik+l-^^P 


T 

£ P k 


-1 


(HO) 
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The vector operators also obey the rule (Ref.1:985): 

r A A 5 Z 

6tOz) = — . z + A • ^ (111) 

T a 

6h (Op 

let A=-— and let z=R~ 1 {r. .. -h( 0 ' )} . Then: 

60" _K 1 K 


P = p" 
k+1 k 


-T 

Ch T (Sp 

/\ 

$ -c 

- r 

A - A 

60 k 2 

r, -h( 0." ) 
—k+1 — k 


6h T (0p 6h(0p 

-- . r -_ 


60 k 


60 k 


T _ 

S P k 


_ 1 


( 112 ) 


Which is equation (54). 

Now substitute equations (25) through (27) into equa¬ 
tions (53) and (54) along with R=o 2 I , where I is a 

n 

P t xP t identity matrix. Write out the vector products as 
summations to give: 


- - 1 ^ ip 

—k+1 — k k+1 — o n Z_/ tt 


sin(p—)sin (p^-( i-l)+§p 


r k+l (i) 


i Ty 

~ sin(^-)cos(^(i-l)+0p 


i T 
s 


(113) 


Utilizing Appendix A 

- - i s TY 

x. = (J )x. - - r-sin 

—k+1 —k TTO 


and some trigonometric identities: 

( fc ) Z> k +i (i)sin( lr (i_1)+ ®k ),p k+i*s 

L L 


(114) 
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which is equation (57). Before doing similarly for equa¬ 


tion (54), move the P R inside the inversion operation on 
the brackets: 


k+1 


-T 

> p k 


0 T - 
« 2 h (O') 
-— R~ 1 

60- 


—k+1 - —^ °k } 


g £ < a P R -1 ^ 9 k> 


60 , 


6 0 ; 


-1 


(115) 


where the fact that P R P R -1 =I was utilized. Now let R * = 

l/o 2 I : 

' n 


k+1 


-T 1 

' P * 


« J h T <e;) 

—|i k+ i-H (S P 




6h(0 R ) 

2 

T 

60 R 


c 


- 1 


(116) 


Substituting for P R from equation (55) and utilizing (25)- 
(27): 

P, 


k+1 


T T -1 T 1 

UP^+GQG 1 ) - c c 1 


n 


i Ty 

-f- ‘ sin( ^> 


,c °s(—(i-l)+ 0 p r k+1 (D- 


i Ty 

s.in(~)cos(^ IL (i-l) 
t t 


+0 k> 


i Tl i Ty 

_§_V( — 

P, Z-P 7T 

t J 


S ) sin 2 (p-)sin 2 (|-^(i-l) + epi 


-1 


( 117 ) 
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The z k (i) defined in equation (63) appears naturally. 

Now substitute this into equation (114), use Appendix A, 

i Y 

and multiply by (—-)(j^)=l . 

s ' 


i TY 

A s 

= *x. - r^r-r 


i Y 


—k+1 — k 7T0 


n 


sln( P7 )( Hr )< rV T 2) z k+i (i) 

t s 


’ sin(—(i-l) + op 


(124) 


W 1 


. , S . -L . . 7T . TT V—* , . „ 

d>x, - (-) •— 2 -*sin(—-).--\ z, ,,(i) 

Y —k tt ' a n p t X S Y k+l v ' 


O tt a 

.sin(^(i-l) + 6p 


(125) 


= (frx,, - 


7(Ft)’ 2 ' 

s 


(i s y)2 T 2 


sin(f-) 

t 


* 2 z k+i (i)sin( r L(i_1)+ ®k ) 


(126) 


2P. 


= 4>x, - 


—k i Ytt 
s 


CNR 


•sin(~)^z k+1 (i)sin(|^(i-l) + ep 


(127) 


The factor CNR was defined in (62) and can be verified by 
substituting equation (15) for and remembering that: 


n'(i) - i n k (i) * «=- - o’ 


(128) 
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appears instead of 


In equation (60) the factor — appears instead of —■ 

s 1 

because it is assumed throughout the paper that the signal 
power is normalized, A-( i s "Y ) 2 = i ; and fluctuations in the 


CNR are modelled only as fluctuations in o 


, a , or 
n ’ n' ’ 


N 0 . The exact same arguments for equation (58) would lead 
to equation (61). 

Now examine the term in equation (127): 


S z k+i (i)sin( l IL(i ~ 1 )+ ®k ) = X> “~ sin (p-) cos (§ JL ( i - 1 ) + e k+1 ) 


sin(|^-(i-l) + 0 p+n^ + 1 (i)sin(|l-( i-l)+ 0 p 


(129) 


= sin (p !! -) sin (0 k - 0 k+1 ) +n k+1 ( i ) sill (|^( i -l) + O k )( 130 ) 


i YP 

-2- Sin( p" ) ^"^k+l^P + n ks 

t 


(131) 


Similarly for equation (118): 


£ z k + l< 1 )c ° s( ?”(i-l) + e') - sl„(l-)cos( 6 k+1 - 0 k ) + n kc 

(132) 
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whore n, is defined in an obvious manner. Since n, 

KC ks 

and n kc are both the sums of weighted gaussian variables, 
they are also gaussian, so it suffices to compute their 


second moments. 


: [" ks ] - E X> 


' +1 (i)sin(|2.(i-l)+0') 


(133) 


2 E [ n k+l (i) sin ( P^ (i_1)+ ®k ) ] 

XX- E „i5[ n i:« (i > sin< p;< i - i )«i;)] 
2>s 


(134) 


(135) 


(136) 


(137) 


Where the identity [x(y)j =E y [ E *|yM] was used, as 

was the fact that n^(i) is zero mean. It is easily shown 
that E [ n kc] = 0 • S ince n^^ and n^ g are zero mean, 

their variances are their second moments. 


[ n ksVs] ■ E SEn- 

_ i J 


(k+1) n ^(k+l)sin(p—(i-l)+0^) 


Ojr /v 

■sin(p-(j-l)+0'. 


(138) 
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= E' 


EE E [n:(k+l)n':(k+l)].§. {cos(|^-(i-l) 
i j n L 1 J -I F t 


A /V 


. 2 TT 


A A 


+O k -0'. )-cos(p-(i+j-2)+0 k +O k .)} 


( 139 ) 


= E; 


£ E n [ n -(k+l)n'(k+l)] -i 


(140) 




kk' 


(141) 


The fact that j^n((k+l)nj(k" + l)j =o^5 kk ^ was used since 

n / is white. Similarly, 


E 


E 




= 0 


(142) 


043) 


So, { n kc .} and {n kg } are s P ectrall y white seqr .^s, at- 
ually independent, and gaussian. 

The above derivations, plus letting (i Y) z /2=l , re- 

O 

suit in equations (65) and (66). 
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Appendix C 

Simulation Subroutine 
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SUBROUTINE GENS IM ( Y ,P.YO » PO » YACT >Q,PHIiKMAX , ITER. 

* IYDIM,IPT.CNRDB,SUMP,SUMY,N.KMPT.WKAREA) 

DOUBLE DSEED 
REAL N 

DIMENSION Y(IYDIM,KMAX),P(IYDIM,IYDIM.KMAX),YO(IYDIM) 
DIMENSION PCX IYDIM,IYDIM),Q<IYDIM,IYDIM) 

DIMENSION YACT(IYDIM.KMAX) ,PH I(I YD IM,I YD IM) 

DIMENSION SUMY(IYDIM.KMAX),SUMP(IYDIM,IYDIM.KMAX) 
DIMENSION WKAREA(KMPT),N(KMPT) 

C 

C *##**#**##****##*##*#### tt**#**#*-#****'#*##**'!!'-*###-*'*#*###*#*' 

C 

C 2D LT MARTIN B. MARK, GE0-80D, 22 SEP 80 
C 

C THIS SUBROUTINE IS A GENERALIZED SIMULATOR FOR THE DE- 
C MODULATION OF A SIGNAL ENCODED ON THE PHASE GF A SINE- 
C WAVE. VIRTUALLY ANY PHASE PLANT MODEL CAN BE ACCGMO- 
C DATED. INTEGRATED SAMPLES OF THE RECEIVED WAVEFORM 
C ARE GENERATED AND USED IN AN EXTENDED KALMAN FILTER 
C TO DEMODULATE THE SIGNAL. 

C THE VARIABLES USED HAVE THE FOLLOWING SIGNIFICANCES. - 
C 

C Y ( I , K > 

C P(I.J.K) 

C YO(I) 

C PO(I.J) 

C YACT(I.K) 

C Q ( I , J > 

C PHI(I,J) 

C KMAX 

C ITER 

C 

C IYDIM 

C IPT 

C 

C CNRDB 

C KMPT 

C 

C SUMP,SUMY, 

C N.WKAREA 

C 

C IMSL LIBRARY ROUTINES GGNML AND LINV3F ARE USED TO 
C GENERATE GAUSSIAN NOISE SAMPLES AND INVERT A MATRIX, 

C RESPECTIVELY. SUBROUTINES NEWP, NEWY, AND SUMM ARE 
C ALSO UTILIZED AND ARE ATTACHED AT THE END OF GENSIM. 

C THE PHASE ANGLE IS ALWAYS THE LAST ELEMENT OF THE 
C STATE VECTOR Y. ALL OTHER ELEMENTS CF THE STATE VEC- 
C TOR ARE ALSO COMPUTED AND RETURNED, BUT THEIR PHY- 
C SICAL SIGNIFICANCES ARE ARBITRARY. 

C 

C**»# ***************************************************** ********* 
INITIALIZE CONSTANTS 


K-TH STATE VECTOR ESTIMATE 

K-TH VARIANCE MATRIX 

INITIAL STATE VECTOR 

INITIAL VARIANCE MATRIX 

ACTUAL STATE VECTOR SEQUENCE 

PLANT NOISE VARIANCE MATRIX 

PLANT TRANSITION MATRIX 

NUMBER OF SEQUENTIAL ESTIMATES OF Y 

NUMBER OF ITERATIONS Y IS AVERAGED 

OVER 

DIMENSIONALITY OF STATE VECTOR 
SAMPLES OF INPUT WAVEFORM PER 
PERIOD 

CARRIER-TO-NOISE RATIO, DB 
NUMBER OF NOISE SAMPLES TO BE GEN¬ 
ERATED BY GGNML 

SCRATCHPAD MATRICES 
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KWD=2*IYDIM 
PI = 3 .141532G54 
PT=FLOAT(IPT) 

CNR =10.**< CNRDB/10.) 

ALPHA=SGRT ( 2 . ) *CNR*PT*-SIN(PI/PT)/PI 
D3EED=19732453GO.DOC 

CLEAR SUMP AND SUMY MATRICES. 

SET FIRST ESTIMATE OF Y TO YO AND 
FIRST VARIANCE MATRIX TO PO 

DO 120 1=1iIYDIM 

DO 110 J=l,IYDIM 
DO 100 K = 1r KMAX 
SUMP <I,J,K)=0.0 
SUMY(J.K)= 0.0 
100 CONTINUE 

P ( I , J , 1 ) =P0< I , J) 

110 CONTINUE 
Y ( I i 1 ) = Y 0 < I ) 

120 CONTINUE 

AVERAGE OVER "ITER" REALIZATIONS OF THE NOISE 
PROCESS BY INCREMENTING DSEED ON EACH PASS. 
NOISE SAMPLES ARE GENERATED BY SUBROUTINE 
GGNML. 

DO 500 1 = 1 ,ITER 
CALL GGNML(DSEED*KMPT 

COMPUTE SEGUENCES OF Y AND P FOR ALL K 


KMAX=KMAX-1 
DO 400 K = 1, KMAX 
KK=IPT*<K-1) 

COMPUTE A-PRIORI UPDATE OF STATE VECTOR 

ARGSUP=0.0 

DO 210 IND=1 ,IYDIM 

ARGSUP = PH I(IYDIM,IND)*Y(IND,K)+ARGSUP 
210 CONTINUE 

SUMS AND SUMC ARE THE CORRELATIONS BETWEEN Z 
AND THE SINE AND COSINE, RESPECTIVELY, OF THE 
A-PRIORI ESTIMATE OF THE PHASE. 


SUMC=0.0 
SUMS=0.0 
DO 300 J=1,IPT 

N(KK+J)=N(KK+J)/SORT<CNR*PT> 

CON=2.*PI*FLCAT(J-l)/PT 
ARG=CON+ARGSUP 

Z=<ALPHA/(CNR*FT)>*COS(CCN+YACT<I YDIM,K)) + 
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* N (KK +J) 

SUMS=Z#SIN(ARG)+SUMS 

SUMC=Z*COS(ARG)+SUMC 

300 CONTINUE 

C 

C UPDATE Y AND P USING SIMULATED DATA AND PLANT 
MODEL. 

CALL NEWP(P,PO,PHI,Q,IYD IM,SUMC,ALPHA,K,KMAX , 

* WKAREA t KUiD . N) 

CALL NEWY(Y,PHI,P,IYDIM-SUMS,ALPHA , K,KMAX » PI) 
CALL SUMM(SUMY,SUMP,Y,P,IYDIM,K ,KMAX) 

400 CONTINUE 

INCREMENT DSEED. 


KMAX=KMAX +1 
DSEED=DSEED+100.D00 
500 CONTINUE 

COMPUTE AUERAGE Y AND P FROM SUMMATION OF ALL 
"ITER" REALIZATION OF NOISE PROCESS. 

DO 620 K = 1 i KMAX 

DO 610 J=1 ( IYDIM 

DO 600 1 = 1 ,IYDIM 

P ( I » J » K )=SUMP(I,J.K)/FLOAT(ITER) 

GOO CONTINUE 

Y(J r K)=SUMY(J » K)/FLOAT <ITER) 

610 CONTINUE 
620 CONTINUE 


RETURN 

END 
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subroutine: newp< P, po, PHI , q, iydim,sumo,ALPHA,K,KMAX, 
-*WKARtA,KWD,S) 

DIMENSION P( IYDIM, I','DIM,KMAX) ,PHI ( IYDIM, IYDIM) 
DIMENSION 0(IYDIM,IYDIM) ,PO(I YDIM,I YD IM) 

DIMENSION WKAREA(KWD) 

C 

C SUBROUTINE NEWP COMPUTES THE NEXT P GIVEN THE 
C LAST P, SUMO, AND THE PLANT MODEL. SUBROUTINE 
C LINV3F IS USED TO INVERT MATRICES WHEN IYDIM 
C IS NOT EQUAL TO ONE. 

C 

DO 400 1=1,IYDIM 
DO 300 J =1,IYDIM 
DUMMYP=0.0 

DO 200 M=l,IYDIM 
DO 100 N=1,IYDIM 

DUMMYP = PH I<I,N)*P(N,M,K)*PHI<J,M)+DUMMYP 
100 CONTINUE 

200 CONTINUE 

PO(I,J)= DUMMYP + Q(I,J ) 

300 CONTINUE 
400 CONTINUE 

IF(IYDIM.EG.1) GOTO 720 
D1 = 1 . 

CALL LINV3F(P0,S,1,IYDIM,IYDIM,D1,D2,WKAREA,IER) 

IF<IER.EQ.130)GOTO 700 

500 P0(IYDIM,IYDIM)=P0<IYDIM,I YD IM)+ALPHA*SUMC 

CALL LIN V 3 F(P 0,S,1,IYDIM,IYD IM,D1,D2,WKAREA,IER) 

IF(IER.EQ.130)GOTO 710 
510 DO G10 1=j,IYDIM 

DO GOO J = 1,1 YD IM 
P(I,J,K + l)=PC(I,J ) 

GOO CONTINUE 

G10 CONTINUE 


700 


C 

710 


C 

C 

rj 


720 


RETURN 

PRINT *, "SINGULAR MATRIX ENCOUTERED IN NEWP" 
PRINT *, -LED BY FIRST OCCURENCE OF LINU3F" 

GOTO 500 

PRINT *, "SINGULAR MATRIX ENCOUNTERED IN NEWP" 
PRINT *, "CALLED BY SECOND OCCURENCE OF LINU3F" 
GOTO 510 


P<1,1,K+1)=1./((1./P0(1,1)>+ALPHA*SUMC) 


RETURN 

END 


uu UUUUCJ UUCJ 


SUBROUTINE NEWY ( Y , PH I , P , I YD IM , SUMS , ALPHA , i< , KMAX , PI ) 
DIMENSION Y(I YDIM,KMAX),PHI(I YDIM,I YD IM) 

DIMENSION P(IYDIM.I YDIM,KMAX) 

C 

C SUBROUTINE NEWY COMPUTES THE NEXT Y GIVEN THE 
C LAST Y, THE NEW P, SUMS, AND THE PLANT MODEL. THE 
C PHASE, Y(I YD IM,K+1), IS RETURNED M0DUL0-2*PI 
C 

DO 200 I=1,1YDIM 
DUMMYY=0.0 

DO 100 J=1,IYDIM 
DUMMYY= PHI( I , J )*Y(J,K)+DUMMYY 
100 CONTINUE 

Y<I,K + l)=DUMMYY-P <I,IYDIM,K+l>*ALPHA*SUMS 
200 CONTINUE 

Y(IYDIM,K+l)=AMOD(Y(IYD IM,K+1),PI) 


RETURN 

END 


SUBROUT INE SUMM(SUMY,SUMP,Y,P,I YDIM,K,KMAX) 
DIMENSION SUMY(I YD IM,KMAX) ,SUMP(IYDIM,I YD IM,KMAX) 
DIMENSION Y<IYDIM,KMAX),P<IYDIM,IYDIM,KMAX) 

SUBROUTINE SUMM ACCUMULATES A SUMMATION OF P 
AND Y FUR ALL "ITER" REALIZATIONS OF THE NOISE 
PROCESS. 

100 DO 300 1=1,IYDIM 

DO 200 J=1,IYDIM 

SUMP(I,J,K)=P(I,J,K)+SUMP <I,J,K) 

200 CONTINUE 

SUMY(I,K)= Y(I,K)+SUMY <I,K) 

300 CONTINUE 

IF(K.LT.KMAX) RETURN 
IF(K.GT.KMAX) GOTO 400 
K = K + 1 
GOTO 100 

400 K =-KMAX 


RETURN 

END 
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